#define  _CRT_SECURE_NO_WARNINGS
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<stdbool.h>
#include<ctype.h>
double length;
double TL, TR;
double den, c, k;
double s;
double dt;
int n, i, tstep, maxtstep;
int tr;
double dx;
double T0[11], T[11];
double ae0[11], aw0[11], ap0[11], ap1[11], b[11];
double ae1[11], aw1[11], x[11];
double a[11], b0[11], c0[11], d[11], p[11], q[11];
void TDMA()
{

}
int main(void)
{
	FILE* output;
	tr = 2;
	length = 3.0;
	TL = 3.0; TR = 5.0;
	den = 100.0; c = 1000.0; k = 10.0;
	s = 10.0;
	dt = 1000.0;
	n = 5.0;
	dx = length / n;
	maxtstep = 2000;
	/*if (tr == 0)
	{
		for (i = 2; i < n; i++)
		{
			ae0[i] = k / dx;
			aw0[i] = k / dx;
			ap0[i] = den * c * dx / dt - ae0[i] - aw0[i];
			ap1[i] = ap0[i] + aw0[i] + ae0[i];
			b[i] = s * dx;
		}
		i = 1;
		ae0[i] = k / dx;
		aw0[i] = k / (dx / 2);
		ap0[i] = den * c * dx / dt - ae0[i] - aw0[i];
		ap1[i] = ap0[i] + aw0[i] + ae0[i];
		b[i] = s * dx;
		i = n;
		ae0[i] = k / (dx / 2);
		aw0[i] = k / dx;
		ap0[i] = den * c * dx / dt - ae0[i] - aw0[i];
		ap1[i] = ap0[i] + aw0[i] + ae0[i];
		b[i] = s * dx;
		for (i = 1; i <= n; i++)
		{
			T0[i] = 3;
		}
		T0[0] = TL;
		T0[n + 1] = TR;
		output = fopen("T.txt", "w");
		for (tstep = 1; tstep <= maxtstep; tstep++)
		{
			fprintf(output, "%d  ", tstep);
			for (i = 1; i <= n; i++)
			{
				T[i] = (ae0[i] * T0[i + 1] + aw0[i] * T0[i - 1] + ap0[i] * T0[i] + b[i]) / ap1[i];
				fprintf(output, "%e  ", T[i]);
			}
			fprintf(output, "\n");
			for (i = 1; i <= n; i++)
			{
				T0[i] = T[i];
			}
		}
		fclose(output);
	}
	if (tr == 1)
	{
		for (i = 2; i < n; i++)
		{
			ae1[i] = k / dx;
			aw1[i] = k / dx;
			ap0[i] = den * c * dx / dt;
			ap1[i] = ap0[i] + aw1[i] + ae1[i];
			b[i] = s * dx;
		}
		i = 1;
		ae1[i] = k / dx;
		aw1[i] = k / (dx / 2);
		ap0[i] = den * c * dx / dt;
		ap1[i] = ap0[i] + aw1[i] + ae1[i];
		b[i] = s * dx;
		i = n;
		ae1[i] = k / (dx / 2);
		aw1[i] = k / dx;
		ap0[i] = den * c * dx / dt;
		ap1[i] = ap0[i] + aw1[i] + ae1[i];
		b[i] = s * dx;
		for (i = 1; i <= n; i++)
		{
			T0[i] = 3;
		}
		T[0] = TL;
		T[n + 1] = TR;
		output = fopen("T.txt", "w");
		for (tstep = 1; tstep <= maxtstep; tstep++)
		{
			fprintf(output, "%d  ", tstep);
			i = 1;
			a[i] = ap1[i];
			b0[i] = -aw1[i];
			c0[i] = -ae1[i];
			d[i] = aw1[i] * T[i - 1] + ap0[i] * T0[i] + b[i];
			for (i = 2; i < n; i++)
			{
				a[i] = ap1[i];
				b0[i] = -aw1[i];
				c0[i] = -ae1[i];
				d[i] = ap0[i] * T0[i] + b[i];
			}
			i = n;
			a[i] = ap1[i];
			b0[i] = -aw1[i];
			c0[i] = -ae1[i];
			d[i] = ae1[i] * T[i + 1] + ap0[i] * T0[i] + b[i];
			i = 1;
			p[i] = -c0[i] / a[i];
			q[i] = d[i] / a[i];
			for (i = 2; i < n; i++)
			{
				p[i] = -c0[i] / (a[i] + b0[i] * p[i - 1]);
				q[i] = (d[i] - b0[i] * q[i - 1]) / (a[i] + b0[i] * p[i - 1]);
			}
			i = 3;
			x[n] = (d[i] - b0[i] * q[i - 1]) / (a[i] + b0[i] * p[i - 1]);
			for (i = n - 1; i >= 1; i--)
			{
				x[i] = p[i] * x[i + 1] + q[i];
			}
			for (i = 1; i <= n; i++)
			{
				T[i] = x[i];
				fprintf(output, "%e  ", T[i]);
			}
			fprintf(output, "\n");
			for (i = 1; i <= n; i++)
			{
				T0[i] = T[i];
			}
		}
		fclose(output);
	}*/
	if (tr == 2)
	{
		for (i = 2; i < n; i++)
		{
			ae0[i] = k / (2 * dx);
			aw0[i] = k / (2 * dx);
			ae1[i] = k / (2 * dx);
			aw1[i] = k / (2 * dx);
			ap0[i] = den * c * dx / dt- k / (2 * dx)- k / (2 * dx);
			ap1[i] = aw1[i] +ae1[i] +ap0[i] + aw0[i] + ae0[i];
			b[i] = s * dx;
		}
		i = 1;
		ae0[i] = k / (2 * dx);
		aw0[i] = k / (2 * dx / 2);
		ae1[i] = k / (2 * dx);
		aw1[i] = k / (2 * dx / 2);
		ap0[i] = den * c * dx / dt - k / (2 * dx) - k / (2 * dx/2);
		ap1[i] = aw1[i] + ae1[i] + ap0[i] + aw0[i] + ae0[i];
		b[i] = s * dx;
		i = n;
		ae0[i] = k / (2 * dx / 2);
		aw0[i] = k/(2 * dx);
		ae1[i] = k / (2 * dx / 2);
		aw1[i] = k/(2 * dx);
		ap0[i] = den * c * dx / dt - k / (2 * dx/2) - k / (2 * dx);
		ap1[i] = aw1[i] + ae1[i] + ap0[i] + aw0[i] + ae0[i];
		b[i] = s * dx;
		for (i = 1; i <= n; i++)
		{
			T0[i] = 3;
		}
		T[0] = TL;
		T[n + 1] = TR;
		T0[0] = TL;
		T0[n + 1] = TR;
		output = fopen("T.txt", "w");
		for (tstep = 1; tstep <= maxtstep; tstep++)
		{
			fprintf(output, "%d  ", tstep);
			i = 1;
			a[i] = ap1[i];
			b0[i] = -aw1[i];
			c0[i] = -ae1[i];
			d[i] = aw1[i] * T[i - 1] + ap0[i] * T0[i] + b[i] + ae0[i] * T0[i + 1] + aw0[i] * T0[i - 1];
			for (i = 2; i < n; i++)
			{
				a[i] = ap1[i];
				b0[i] = -aw1[i];
				c0[i] = -ae1[i];
				d[i] = ap0[i] * T0[i] + b[i]+ ae0[i] * T0[i + 1] + aw0[i] * T0[i - 1];
			}
			i = n;
			a[i] = ap1[i];
			b0[i] = -aw1[i];
			c0[i] = -ae1[i];
			d[i] = ae1[i] * T[i + 1] + ap0[i] * T0[i] + b[i]+ ae0[i] * T0[i + 1] + aw0[i] * T0[i - 1];
			i = 1;
			p[i] = -c0[i] / a[i];
			q[i] = d[i] / a[i];
			for (i = 2; i < n; i++)
			{
				p[i] = -c0[i] / (a[i] + b0[i] * p[i - 1]);
				q[i] = (d[i] - b0[i] * q[i - 1]) / (a[i] + b0[i] * p[i - 1]);
			}
			i = n;
			x[n] = (d[i] - b0[i] * q[i - 1]) / (a[i] + b0[i] * p[i - 1]);
			for (i = n - 1; i >= 1; i--)
			{
				x[i] = p[i] * x[i + 1] + q[i];
			}
			for (i = 1; i <= n; i++)
			{
				T[i] = x[i];
				fprintf(output, "%e  ", T[i]);
			}
			fprintf(output, "\n");
			for (i = 1; i <= n; i++)
			{
				T0[i] = T[i];
			}
		}
		fclose(output);
	}
}
